Fictitious excitations in the classical Heisenberg antiferromagnet on the kagome lattice 
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Using an advanced Monte Carlo algorithm and high precision spin dynamics simulation, we inves- 
tigated the dynamical behavior of the Heisenberg antiferromagnet on a kagome lattice at extremely 
low temperatures. We demonstrate that the detection and correct identification of propagating 
excitations depends crucially on the choice of coordinates and we show how modes are displaced in 
Fourier space if a single, global Cartesian coordinate system is chosen. 



I. INTRODUCTION 

Although frustrated magnetic systems have enjoyed 
the attention of the scientific community for many years, 
they remain fascinating research subjects and their rich 
behavior is far from being completely understood. If only 
classical interactions are considered, it is usually an easy 
task to find the ground state, which is often degenerate. 
But once the systems are thermally excited, intriguing 
effects such as vortices^ or magnetic monopoles^ can be 
observed. 

Computational physics provides important tools for 
probing the behavior of these magnets. Furthermore, it 
bridges the gap between theory and experiment by test- 
ing theoretical predictions which are usually stated for 
idealized systems that do not exist in their pure form 
in nature and are, therefore, not accessible by experi- 
ments. On the other hand, more realistic models are 
often too complex for rigorous analysis but are to some 
extent treatable in computer simulations, which are be- 
coming increasingly powerful due to improved algorithms 
and faster machines. Spin dynamics methods are valu- 
able in this respect because they allow the determination 
of the dynamic structure factor which is closely related 
to results of inelastic neutron scattering experiments. 

One system of great interest is the kagome Heisen- 
berg antiferromagnet (KHAFM). Studies of quantum an- 
tiferromagnets helped pique interest in corresponding 
classical systems, see Ref. and both experiments 

and numerical studies of quantum systems continue to 
this day^i^. Twenty years ago an extensive theoretical 
description^ii was developed for the classical Heisen- 
berg model, and a number of numerical studies^i"— have 
been performed since. Predictions for dynamical behav- 
ior were, in part, confirmed by recently published results 
of spin dynamics simulationsii. However, a branch of 
unprcdictcd modes was found and its origin has not yet 
been clarified. This makes a closer look at the subject 
worthwhile. 

The paper is organized as follows: After this introduc- 
tion we will introduce the model and review some es- 
tabhshed facts about its behavior. In section HI, the 
methods used in our work for simulation and analy- 
sis are presented. Section IV is devoted to a presen- 
tation and discussion of the results of the Monte Carlo 
and spin dynamic simulations, followed by concluding rc- 
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FIG. 1: The kagome lattice and the two lattice vectors ai and 
SL2. Spins are placed at the line intersections. The [10] and 
[11] lattice directions are shown by (unit) vectors eio and en. 



marks in Section V. The paper concludes with two ap- 
pendices in which we discuss some calculations in greater 
detail and introduce a modification of the Wang-Landau 
algorithmi2r2i. 



II. BACKGROUND 

In the classical Heisenberg model, three dimensional- 
vectors s of unit length interact via the Hamiltonian 



n 



(ij) 



(1) 



with a negative interaction constant J favoring antifer- 
romagnctic spin pairs. The kagome lattice is drawn in 
Fig. [TJ We apply periodic boundary conditions and as- 
sume in the following that the lattice constant |a,i| ~ 1. 
If N is the number of sites, the system contains 2N/3 
triangles. The energy is minimal (i^min = — A^|>^|) if the 
sum of the three spins belonging to each triangle is zero. 
This leaves the system underdetermined (frustrated) and 
the ground state is highly degenerate. 

At low enough, but non-zero, temperatures a common 
spin plane is established, at least in finite systems, in a 
process called 'order by disorder—. This 'coplanar' state 
is stabilized by the entropy of local modes which can be 
constructed on the hexagons of the kagome lattice. Only 
three distinct spin direction occur separated by angles of 
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FIG. 2: (a) Typical configuration at low temperatures. Spins 
fluctuate around three distinct directions which are indi- 
cated by different symbols and colors, (b) Solid lines show 
{weathervane loops which are constructed if spins of two 
'types' are connected. 



2tt/3 which allows the assignment of labels ct; S {1,2, 3} 
to each spin according to the direction towards which 
it is oriented. 

It is then possible to classify configurations {s} ac- 
cording to their labels {a}. The ground state condition 
requires that each triangle has spins with three differ- 
ent labels; hence, valid {a} arc also ground state con- 
figurations of a three-state Potts model. However, while 
the probability distribution over all {a} is uniform in 
the case of the Potts model, there is a selection in the 
KHAFM due to cntropic effects. This can be understood 
if the concept of weathervane loops is employed. These 
loops form when two label values are chosen and spins 
possessing these values are connected along the bonds 
of the lattice (Fig. [Hb]). Because there are three ways 
of choosing two out of the three possible values, three 
overlapping loop structures coexist and each spin belongs 
to two loops. However, within a single structure, loops 
never touch but are separated by spins pointing in the 
third direction. Therefore, spins of a single loop can ro- 
tate collectively around this direction while causing only 
small changes to spin-spin angles and, thus, to energy. 
However, unlike in the case of a true weather vane com- 
plete rotations rarely occur. Instead, relatively large fluc- 
tuations in these collective degrees of freedom occur and 
account for the difference between the ground state of 
the Potts model and the distribution over different {a} 




FIG. 3: (a) The VSxVS structure; (b) Weather vane loops for 
the VSxVl structure are shown by solid lines. 



in the KHAFM. Because each loop provides one degree 
of freedom for fluctuations, a structure with many loops 
has a higher entropy than a structure with few loops. In 
a counter effect, however, entropy decreases if the loop 
number becomes very large, and the maximum loop num- 
ber does not occur in large enough systems, because it is 
realized in only six distinct configurations {cr} possessing 
the so-called -\/3x "v/S structure (Fig. In the -\/3x "v/S 
state each hexagon is a weathervane loop and the loop 
structure is periodic in three directions with the six wave 
vectors k^^^g = ±|7rai, ±|7ra2, ±|7r(ai - aa). 

In an analytical study,— Harris et al. examined the 
KHAFM with and without- second and third-nearest- 
neighbor interactions and predicted a branch of soft 
modes with frequencies approaching zero as the temper- 
ature decreases as well as a two-fold degenerate acoustic 
branch. In the case of only nearest neighbor interactions 
the frequencies of the acoustic branch were given by 



a;(k) = |J| J2(sin2qi -fsin^ga +sin2(qi -92)), (2) 



where qi = k,^^ q2 = {kx — V3ky)/2, and k is the wave 
vector. 
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FIG. 4: With decreasing temperature spins tend towards three 
basic directions and labels, a, can be assigned to each spin ac- 
cordingly. Fluctuations around these directions are described 
in the alternative coordinates u, v given in Eq. [7) 



III. METHODS 

A. Alternative Coordinates 

Because no anisotropy and no external field are con- 
sidered, the Hamiltonian lias spherical symmetry in spin 
space and the system can be freely rotated. In the fol- 
lowing, we will choose the orientation in such a way that 
the spin plane coincides with the xy-planc. and wc name 
the angle between the average direction of the a = 1 
spins and the y-axis (j>. Alternative in-plane coordinates 
s", s"" , as used by Harris et al.— are now given by 



I cosi/'i sinipi 
$1 = — sin cos 
sf / \ 1 




(3) 



where = + 27r/3(cri — 1) is the angle enclosed by 
the j/-axis and the basic spin direction which is denoted 
by Gi. The new it-axis is perpendicular and the w-axis is 
parallel to this direction (Fig. 0]). 



B. Monte Carlo simulations 

We used the simulated temperingi^ technique which 
allowed us to cover a large temperature range. It pro- 
vides a generalized ensemble that is a composition of Nt 
canonical ensembles with a priori defined temperatures 
Tj. The probability to find the system in a configuration 
{s} at temperature Tj is 



..({s}) = — expf-^ 



(4) 



where Wi is a weight factor and the normalization con- 
stant ZgT is the partition sum of the entire ensemble 



^ST 



Nt 



(5) 



i=l 



which is here written as a weighted sum of canonical par- 
tition sums Z . The acceptance probability for a Monte 
Carlo step from configuration {s} to configuration {s'} 
and temperature Ti to temperature Tj is given by 



/^-({s},{s'}) = min(l,^exp 



H({s}) H{{s'}) 



(6) 

It follows from Eq. ^ that conformational updates 
can be performed just like in a Metropolis simulation as 
long as the temperature remains unchanged. Usually, 
in order to move in temperature space, separate Monte 
Carlo steps which do not alter the configuration of the 
system are employed. To ensure that all temperatures 
are sampled with equal probability, the weights have to 
be chosen such that Wi oc Z{Ti)^^. Since the partition 
function is not known at the beginning, we performed 
preliminary simulations to determine Wi using a modifi- 
cation of the Wang-Landau algorithm that is described 
in Appendix A. For our simulation, we used Nt ~ 10000 
different temperatures distributed equidistantly on a log- 
arithmic temperature scale: —6 < logiQ{kBTi/\J\) < 3. 
To update configurations at constant temperature, we 
mainly used the heat-batbi^ method; in doing so, we 
avoided the problem of too low or too high acceptance 
rates as result of improperly chosen step length. Details 
of the heat-bath technique for a Heisenberg spin system 
have been published by Miyatake et alJ^. 

Obtaining thermodynamic quantities from simulated 
tempering simulations is relatively simple. Although the 
rcweighted histogram techniquei^ could be applied, it is 
not required here. It follows from Eq. (0]) that the dis- 
tribution P{E) can be written in terms of a density of 
states g{E): 



Nt 

P{E) cx g{E) exp 



1=1 



(7) 



Once the weight factors are chosen it is easy to calculate 
the sum on the right hand side, and the density of states 
can be estimated via a histogram over E. 

As described above, the coplanar state is highly degen- 
erate. Different configurations {ct} exist in great num- 
ber and, although requiring no increase in energy, transi- 
tions between them are inhibited by cntropic bottlenecks. 
Hence, autocorrelation times are large in the coplanar 
state and ergodicity might even be broken at very low 
temperatures. This is partially corrected for by use of 
a generalized ensemble method, such as simulated tem- 
pering in our case, but further improvement is possible 
if jumps between different configurations a are enabled 
at low T. To do this, we introduce a Monte Carlo step 
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FIG. 5: Spins of a weather vane loop can be updated collec- 
tively in a reflection. The boundary spins (red) remain un- 
changed. The mirror plane (rectangle) contains the average 
direction of the boundary spins and is perpendicular to the 
(local) spin plane. 



that is designed to "flip" all spins along a weathervane 
loop. Here, "flipping" means exchanging spin directions 
(T that already occur in the loop, i.e., if a weather vane 
loop is composed of spins with numbers < Moop and 
CT/. = 1 or CT;. = 2 for all spins in the loop, then in the 
modified loop it is still cr^ = 1 or cr[ = 2 for all i but at 
exchanged positions cr^ ^ cr,- (FiglS|). A local reference 
frame is established with the use of the loop spins and 
the neighboring (boundary) spins hi, for which in this 
example tTf,. =3. 

The first step is the identification of a weather vane 
loop. While this can be done on a global scale by deter- 
mination of the configuration {u}, we used a more flexible 
procedure that, in principle, can result in accepted up- 
dates even if the system is in a disordered state. In the 
beginning, one spin is randomly chosen as the first spin of 
the loop sjj and one of its neighbors is - again randomly 
- drawn as the second . Their common neighbor is rec- 
ognized as a boundary spin s^^ which is not, and cannot 
be, part of the loop. To continue the construction of the 
loop, the following procedure is repeated until the loop is 
closed; Assuming that n spins are part of the incomplete 
loop, then one neighbor of s;^ is part of the loop (sj^^ J 
and one neighbor is the boundary spin Sbr^-i- The next 
spin (s;^^j) must, thus, be chosen from the two remaining 
neighbors s^i and s^a, which arc neighbors themselves. 
Given the nature of a weather vane loop, this spin has to 
be oriented in approximately the same direction as s;^ j, 
so In+i = mi if s;„_j • s,„j > S;„_j • s,„2 and /„+i = m2 
otherwise. The spin that is not chosen is consequently 
the boundary spin s;,^ . If during this procedure a spin is 
included that is already a boundary spin or part of the 
loop but not , the procedure has failed and the Monte 
Carlo move is rejected. 

There are two intuitive possibilities to update the loop. 
The first one is a rotation of the spins of the loop in 
spin-space by an angle of tt around the direction given 
by the boundary spins. This, however, also alters the 
orientation of the spins with respect to the spin plane. A 
spin that is pointing above the plane will afterwards point 
below it and vice versa. This is not desirable because each 
spin also belongs to a second loop, and the out-of plane 



component of the spin coordinates is mainly determined 
by the collective excitation of each loop. Hence, a switch 
in this component distorts the second loop and increases 
the energy. One has to choose the alternative, which can 
be described BjS 0, refiection on a plane spanned by the 
direction of the boundary spins and the normal vector 
on the spin plane leaving the out-of plane component 
untouched. 

If a closed loop could be constructed, then two vectors 
V;,, v; spanning the local spin plane could easily be found. 
The first one is the sum of all boundary spins 

^loop 

V6 = ^ s,, (8) 

1=1 

and a roughly perpendicular second vector is created by 
alternating the addition and subtraction of the loop spins 

A^ioop/2 

v; = XI -Si2.-i- (9) 
1=1 

A normal vector of the plane on which the spins are re- 
flected is then given by 

u = Vb X (vb X Vi), (10) 

and the updated spins arc 

s;^=s,,-2u^. (11) 
|u| 

After each loop flip we must test if the above procedure 
would also identify the flipped loop. If this is not the case, 
the inverse move is impossible and the update has to be 
rejected. The spin plane and the mirror plane will always 
be identified correctly in the inverse update. Finally, the 
change in energy has to be determined, and the move is 
accepted with the probability given in Eq.Q. 

It turns out that in the coplanar state, the acceptance 
rates for this update are independent of temperature and 
decrease exponentially with increasing loop length. This 
suggests that the average change in energy is positive and 
proportional to the loop length. 

C. Spin dynamics simulations 

Employing the described Monte Carlo algorithm, we 
extracted configuration {s} from canonical ensembles at 
temperatures Ti . The coupled equations of motion 

s,(t) = JH,{t) X s,(i) (12) 

were then integrated over a time interval oit — 2000| J|^^ 
using a 4th order Adams-Bashforth-predictor-Adams- 
Moulton-corrector methocU^ with a time step of At = 
10~"'| J|~^. Here, the local magnetic field Hi is the sum 
of the spins adjacent to Sj. The resulting trajectory is 
subsequently used to calculate space displaced and time 
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FIG. 6: Probability pl that a constructed loop is of length L, 
and probability p'^'^ that an attempted flip of a loop of length L 
is accepted. Data in the lower right represent loops that span 
the entire system which for this case (24? x 3 spins) requires 
L > 24. (Note that the number of spins m a loop Nioop ts 
twice its length L in units of \aLi\.) 



displaced correlation functions, e.g., correlations of the 
spin's x-component: 



'(0) 



(13) 



where r' runs over all lattice sites. Finally, a space-time 
Fourier transform results in the dynamic structure factor 



(14) 



(5'^,S'^ similarly), which in turn is closely related to 
outcomes of inelastic neutron scattering experiments, if 
global coordinates (x, y, z) are used. For the KHAFM, 
the more genuine description, however, is given by the al- 
ternative coordinates (w, u, z) and we also calculated the 
respective structure factors (5", S*"). 

Spin dynamics methods have proven to be valuable 
for the study of diverse dynamic excitations, e.g., spin- 
waves^S, vorticesi, and solitons^i. 



IV. RESULTS AND DISCUSSION 

A. Static Properties 

As a first result we obtained the density of states g{E) 
and calculated the specific heat (Fig. [7]). At low T on each 
of the iV/3 hexagons a local soft mode can be constructed. 
These modes experience anharmonic potentials of fourth 
order with a thermal energy of i/cBT. All other de erees of 
freedom contain ^k^T , leading to a total average energy 



C 




-0.2 



-0.4 



-0.6 



N\J\ 



E/N\J\ + l / 







10-6 IQ-'^ 10"* 10-3 10-2 10"^ 1 
kBT/\J\ 



10 



of {E)/N\J\ 



11 



fceT — 1 at small T. As in previous 



FIG. 7: Specific heat C/N , mean energy {E)/N, and density 
of states g{E) for a system containing N = 6912 (— 48'^ x 3^ 
spins. Statistical errors do not exceed the line width. 



studies^iSiiS this theoretical prediction characterizing the 
coplanar state is fulfilled for T < 3 • IQ-^] J|/fcB. 

At intermediate temperatures, C/Nk^ ~ 1 and low 
E indicate that only short range correlations exist and 
that spherical symmetry is not yet broken in favor of a 
common spin plane. 

Note the linear shape of the density of states in the 
double-logarithmic plot in the inset of Fig. [7l Although 
of minor relevance here, this behavior is a general feature 
of classical systems corresponding to a finite limit of the 
specific heat for T — > 0. This can be shown easily using 
the following ansatz for the density of states: 

\og^^j\og{{E-Eo)/\J\), (15) 

where Eq is the ground-state energy Eq ~ — A^|'/| and G 
is a normalization constant. One finds that the specific 
heat C = ^k-Q and the mean energy {E) = Eq + ^k^T . 
A fit of the data in Fig. [7] for \og^QkBT/\J\ < -5.5 
gives 7 = 6335.9(2) which is close to the value 
N ■ 11/12 = 6336 which we expect for T ^ 0. 

From simulation and theory^^i^i^ it is known that the 
probability distribution over different configurations {a} 
is non- uniform and biased toward the -\/3x-\/3 state. 
It has been under discussion whether this distribution 
changes with temperature and if long-range order is es- 
tablished for r — )• 0. To investigate this possibility and 
to test the algorithm, we produced thousands of confor- 
mations of a 432 spin system at different temperatures, 
cooled them below T = 10^'''| Jj/fce, and performed sim- 
ulated tempering simulations without loop flip updates 
at IQ-'' < kBT/\J\ < 10"^. Due to the low tempera- 
ture, ill more than half of the simulations no loop flips 
occurred spontaneously, which enabled a comparison of 
average logarithmic temperatures of single conflgurations 
{cr}. If certain conflgurations are favored at low temper- 
atures, then this would result in values closer to the lower 
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FIG. 8: If the simulation is constrained to logjQ(fcBT/| J|) < 
—5 the time-averaged logarithmic temperatures are indepen- 
dent of the original temperatures Tinit- The broken line is a 
guide for the eye. 



bound while configurations taken at higher T would con- 
sequently show higher average logarithmic temperatures. 
The results are shown in Fig. [H and within the precision 
achieved no evidence for such a trend can be found. The 
exact values have no physical meaning; with perfectly 
tuned weights Wi the distribution over the logarithmic 
temperature would be flat and the averages would con- 
verge to 5.5. 

We conclude that for log]^g(A;Br/| J| < —3 the ensem- 
ble of loop structures is practically independent of tem- 
perature. Consequently, no further fundamental order is 
established and changes in energy, sub-lattice magneti- 
zation, or other order parameters'^ are based solely on 
decreasing excitations on a fixed set of ground-state con- 
figurations. 

It should be noted, however, that Henley^ recently 
studied effective Hamiltonians and suggested that long 
range order may exist for T — )■ at length scales that 
exceed the size of the systems we could investigate and 
which could, therefore, be neither confirmed nor excluded 
by our simulations. 

It is instructive to investigate the static correlations 
between labels cr, which we define as: 



(16) 



This function bears no noticeable temperature depen- 
dence as long as the system is in the coplanar state. 
This confirms once more that the distribution over the 
loop structures does not change with T. Yet, Fr is best 
measured at lower temperatures where {cr} can be de- 
termined most easily. As previous studies observed, we 
find that the systems show -s/S x -\/3 correlations which 
decrease with increasing distance. The data are well de- 
scribed by a power law: 



where 



rf^^=^os(rk^,^). 
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FIG. 9: Static correlations Fr for a 48^ X 3 lattice follow a 
decaying p/^^^ pattern if the system occupies the coplanar 
state. The solid lines show a power law oc |r|~"'^''^^® with in- 
cluded periodic boundary conditions. (In [11] direction Fr has 
a period o/ • 48 ^ 83.1 ) 



To estimate the exponent k, we fitted a power-law func- 
tion G with included periodic boundary conditions. 



G(r) = Go ^ 5] |r - zLai - jLaap'^ 



(19) 



to Fr/F^'^^ and obtained k w 1. Some of the relevant 
data and the fit for a system with L = 48 are shown in 
Fig. [91 



B. Dynamic properties 

In our investigation, we consider solely wave vectors 
parallel and perpendicular to the bonds of the system, 
which simplifies Eq. © to 



uj{lew) = |J|V3 - 2 cos Z - cos2;, (20) 
with the unit vector eio = ai, a2, ai — a2 and 



c^(/en) = 2|J| 



sm 



(21) 



with eii 



2ai— a2 



-ai+2a2 



lai+aal' |2ai-a2|' |-ai+2a2| 



(Fig.P). 



All data presented in this section are from simula- 
tions of a system containing 48^ x 3 spins at T = 
lO^^I J|/A;b- We produced 2457 independent configura- 
tions and performed spin dynamics simulations of length 
tsD 2000| J|-i with a step length of M = IQ-^I J|-^ 
Cauchy distributions were fitted to the maxima, and the 



FIG. 10: Dynamic structure factor m the (10) direction based 
on the alternative in-plane coordinate u. Maxima agree well 
with the dispersion relation from Eq. f2U\) which is plotted 
as a line in lu) -plane. Due to the low temperature systematic 
deviations (inset) are very small. For the sake of convenience, 
we increased S by a factor of lOT" here and in Figs. and 

[21 



jackknife method was employed to estimate statistical er- 
rors in peak positions. While we will focus on the modes 
of the acoustic branch, we acknowledge that prominent 
signals at a; = give strong evidence for soft modes. 

We also obtained results for smaller lattices and ob- 
served the systematic variation in the dynamic structure 
factor with lattice size. The results for smaller lattices 
were useful for the eventual interpretation of the large 
lattice data but do not themselves show different phe- 
nomena. For this reason, we do not show those data 
here. 

We find that the results for alternative coordinates 
u,v,z agree very well with the predictions. In both direc- 
tions, the soft modes produce prominent peaks at very 
low frequencies. While no trace of acoustic modes can be 
found in S"", prominent maxima exist in (Figs. [TOIlip 
and S*^ at the frequencies given by Eqns. (|20l2ip . 

The in-plane dynamic structure factor in global coor- 
dinates = S"^ -I- S*^, however, shows a completely 
different behavior. In the (10) direction (Fig. I12[a]) arcs 
of maxima intersect, but no excitations of comparable 
intensity occur at the frequencies of the acoustic branch 
according to Eq. ([20]) . In the (11) direction an almost 
dispersionless branch of maxima with high frequencies 
resembles an optical mode and was, in fact, interpreted 
as such by Robert et aliii. 

Additionally, numerous ridges of low intensity form 
washboard-like patterns in S^^ in both directions, but 
these vanish if a pure structure is considered 

(Fig. I12[b]). This washboard pattern is pronounced for 
small lattices but washes out as the lattice size increases. 

To understand this behavior, we need to know the re- 
lationships between the in-plane correlation functions in 



FIG. 11: Dynamic structure factor for the alternative in-plane 
coordinate u in the (11) direction. The behavior is well de- 
scribed by the dispersion relation m Eq. i21p (solid curve). 
Deviations of peak positions are small (inset).) 



global and alternative coordinates. In Appendix B we 
show that 

C^y{t) (xTrC^it). (22) 

It follows that C^y{t) and S'=y{k,uj) depend on the un- 
derlying loop structure {a} even if C" (i) does not. For 
the "v/S X a/3 structure the Fourier transform of r/*^^ is 

r^^^(k)-E^(k-kV3,^), (23) 

where S is Dirac's delta. Using the convolution theorem 
we find 

6 

5-^(k,c.)cx^5"(k-k*^^,a;). (24) 

This means that modes seen in iS" contribute six-fold to 
S"^^ as fictitious modes displaced by k^^^y^ as shown in 

Figs, [m and [T5l The particular choice of k^^^ is ir- 
relevant unless intensities are considered. Due to the pe- 
riodicity of w each of the three directions will cause the 
same apparent branches by shifting the center of the Bril- 
louin zone onto a corner and the opposite corner into 
the center. In the (10) directions this results in acoustic 
branches which follow the same dispersion relation as Eq. 
(PO]) but shifted by ±27r/3. 

In the (11) direction, the apparent optical mode turns 
out to be a part of the acoustic branch and the dispersion 
relation (for the shifted positions) reads 

(I;(;eii) = a;(/eii ± k^^^), 

= y^^-f cos (VSZ). (25) 
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(a) 




FIG. 12: Dynamic structure factor m the (10) direction based 
on global in-plane coordinates. A shift in Fourier space oc- 
curs (see text) and no signals matching the analytical disper- 
sion relation (continuous lines) are observed. In the canonical 
ensemble (a) a washboard-like pattern is caused by the power 
law-like decay of the \/3x-\/3 correlations; (b) In case of long- 
range correlations, e.g. for the pure structure, this 
pattern vanishes. 



The described behavior is induced by the VSx a/3 cor- 
relations which in the general case of free {cr} are present 
but appear to decay in a power law (Fig. ^ . This means 
that we have to apply the convolution theorem a second 
time to describe the shifts in fc-space if an ensemble of 
typical loop structures (not pure -y/Sx'v/S) is considered. 
Since the Fourier transform of the power law is again a 
power law, it follows that in principle each mode can be 
measured at each point in Fourier space with maximal in- 
tensity at the position given by the initial shifts ztkyg^yg. 
This effect causes the washboard pattern by smearing 
the peaks in k-space; however, as mentioned earlier, the 
shifted peaks can only be observed for finite systems. 
We expect no signs of acoustic modes in ^^'^(k, w) in the 
thermodynamic limit. If, on the other hand, -s/Sx-s/S 
correlations are stabilized by second- and third-nearest 
neighbor interactions^ the shifted acoustic branch with- 
out the washboard would be detectable in experiments 




FIG. 13: As in Fia. \lS\. a shift in Fourier space can also be ob- 
served m the (11) direction, and modes from a different part of 
the Brillouin zone manifest at frequencies given by Eq. k25\) , 
drawn as a broken line. Peak positions used for calculating 
the deviations in the inset were determined from simulations 
of the VSxVS structure allowing a higher precision. 



u/\J\ 

0.5 1 1.5 2 2.5 




k, 



-2 2 4 




FIG. 14: Analytical frequencies of the acoustic branch as a 
function of wave vector k. In the (10) direction, the shifts 
by kyg^yg and by — ky^^yg contribute differently, resulting in 
two apparent branches. The gray hexagon marks the edge of 
the first Brilloum zone. 
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|ai-a2| 
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FIG. 15: Analytical frequencies of the acoustic branch as a 
function of wave vector k. In the (11) direction, the shift 



by k 



creates a fictitious optical branch. The shift by 



-k^^^ (not shown) has an identical effect. 



and might actually have been observed allready^. 

Considering the described displacements, we find that 
all observations become consistent. Not only do the po- 
sitions of maxima in 5"(k, w) and 5"^^(k, w) agree nicely 
with theory, the deviations caused by non-zero tempera- 
ture also match. From Fig. [TU] and Fig.[T3]it follows that 
at the corners of the Brillouin zone, acoustic modes have 
frequencies above the analytically calculated values while 
no deviation is noticeable at the midpoints of the zone 
edges according to Fig. [TU] and Fig. [TTJ These two plots 
also show that close to the center, frequencies are slightly 
increased, while they are below theoretical predictions at 
intermediate distance from the center (for the latter, see 
also Fig. 1131). 



orders of magnitude lower than previously reached. To 
do so, we designed a trial move that allows the "flipping" 
of weather vane loops, and we combined the simulated 
tempering and heatbath techniques with a modification 
of the Wang-Landau algorithm. We were able to show 
that, in agreement with theory, minimal temperature- 
dependent selection of loop structures occurs over three 
orders of magnitude in temperature. 

Extensive spin dynamics simulations at k-BTi/\J\ = 
10~^ provided results that agree very well with analyti- 
cal calculations: Using suitable coordinates, we observed 
harmonic modes with the predicted frequencies directly 
within the spin plane and transverse to it. Considering 
global Cartesian coordinates, we showed that signals in 
the dynamic structure factor are shifted due to a/Sx \/3 
correlations. Apparent in-plane excitation at the center 
of the Brillouin zone with non-zero frequency actually 
belong to acoustic modes at the zone edge. Hence, the 
classification of these signals as a branch of optical modes 
by Robert et al.— was premature. 

We argue that due to the decay of the correlations, 
in-plane excitations will not be measurable for large sys- 
tems unless next-nearest neighbor interactions stabilize 
the -\/3x correlations. 



VI. ACKNOWLEDGEMENTS 

We thank Shan-Ho Tsai for helpful discussions. This 
work was funded by NSF grant No. DMR-0810223. 

VII. APPENDIX A 

The Wang-Landau algorithmi^r— is a method that is 
designed to adjust a weight function Wi{Q) over time i 
such that a flat histogram H{Q) can be produced. The 
parameter Q is often the system energy which allows 
the determination of thermodynamic quantities such as 
specific heat and mean energy via the estimation of the 
density of states as a function of energy. Because we 
are interested in producing a flat histogram over a log- 
arithmic temperature range, we will treat Q as a gen- 
eral variable in these considerations. In the simplest 
case, the probability of a configuration X is desired to 
be proportional to w (Q(X)) which means that proposed 
Monte Carlo moves X — > X' are accepted with probabil- 
ity Pacc(X,X') = minil, w{Q{X')) /w{QiX))). In that 
way a series of configurations X.; and values qi — Q(X.i) 
is produced. In order to achieve fiatness we modify the 
weight function after each step: 



V. CONCLUSIONS 

For this work we devised a Monte Carlo algorithm that 
was able to equilibrate the Heisenberg antiferromagnet 
on the kagome lattice at temperatures more than two 



Wi+i{Q) = 



Wi{Q) otherwise 



(26) 



where / is a modification factor which is initially set to 
/ = e (Euler's constant) and gradually reduces to unity. 
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thus causing decreasing changes to w. We propose a mod- 
ification of the original Wang-Landau algorithm which 
is based on the assumption that especially early in the 
iteration, the deviations from the detailed balance crite- 
rion can be reduced by delaying the modification of the 
weight function. It is self-evident that a change in w af- 
fects the simulation's balance most if it is performed in 
the proximity of the current position qi of the system. 
We, therefore, introduce a delay for the modification of d 
time steps which allows the simulation to proceed undis- 
turbed, decorrelate, and be less influenced by alterations 
of w. When a modification at qi-d is eventually made, 
it is necessary to take into account the changes of w at 
that position between time i — d and i. The corrected 
modification formula reads: 



w,+,{Q)=l . (27) 

I, 'Wi{Q) otherwise 

Note that these equations become equal to Eq. [26] if d = 
0, i.e., in the case of no delay. To calculate this formula, 
it is necessary to know qi-d and Wi-d{qi-d), i-c at each 
time i the values qi and Wi{qi) have to be stored for d time 
steps giving a total of 2c? numbers to be held in memory. 
On modern computers, this can be done for large values 



of d; however, we found that a good compromise between 
increased balance and fast convergence is achieved at c? sa 
103,10". 

It turns out that systematic errors are much smaller 
in the beginning of the simulations and that the overall 
convergence is often faster, but in our experience never 
slower, than in the original Wang-Landau method. In 
the past this method was used successfully to calibrate 
one- and two-dimensional weight functions for polymer 
simulations^!^. 

In the present study the weight function w corre- 
sponds to the weight factors Wi, and the quantity Q is 
logj^g fceT/l J|. Modifications of W were performed after 
each temperature update. 



VIII. APPENDIX B 

Consider the product of the x-components of two spins 
as a function of local coordinates s". The three main spin 
directions lie in the xy-plane, while one of which and the 
y-'Axis include an angle (j) (Fig. [4]). Instead of taking the 
average over all possible angles 4>, we average only over 
circular permutations of spin directions cr, which corre- 
sponds to rotations of all spins by angles ±|7r around 
the z-axis. 



If cr,; = (To 



If ^ aj-. 



3 ■ 



cos 



COS 



-TT + COS 



4 4 
. s~. I 3 + COS 2(p + cos(20 -I- -tt) -|- cos(2(/) — -tt) 

o o 



(28) 



-s^.s^. I cos0cos( 



-TTj + COS(pCOS( 
O 



2 2 2 

-tt) + cos{(f) + -tt) cos((/) - -tt) 



1 / 2 2 2 2 4 

-Sr. Sr. cos(20 -I- ^tt) + cos(--7r) + cos(20 - -tt) 4- cos(-7r) + cos(20) + cos(-7r~ 
D • ^ \ 3 3 3 3 3 

1 / 2 2 4 \ 

gSr.Sr, ( cos(--7r) +cos(-7r) -|-cos(-7r) I , 



(29) 



In both cases, the contributions of (j) vanish which 
means that the average is equivalent to an average over 
(j). Furthermore, {Sr.Sr^)^cr} = (^ri*rj){cr} even if is 
constant because an average over all loop structures {cr} 
involves the average over permutations. Following Harris 
et ali^ we exploit the fact that collective harmonic exci- 



tations s" are independent of the spin configuration {a} 
and average over all {cr}: 
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3p _ 1\ ^onA It follows that 

3 

where fri-rj denotes the probability for ai = aj-. (0 ^ ^TrCj. {t) (33) 



fr.-r, = (31) 



and 



with 5 being the Kronecker delta. Inserting the afore- 
mentioned correlation function F, we obtain C^y{t) = C^{t) + C^{t) = — rrC"(t) (34) 



2 



«4)W = isrVr,rr.-r,. (32) 
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